;+
;
; NAME:
;    
;       
; PURPOSE:
;    Find the elements of an array of integers (arr2, which may contain
;    duplicates) which have a match in another array (arr1, which may not
;    contain duplicates).
;
;    No attempt is made to find the matches from arr1.  Use the program
;    match_dup for that purpose.  This program is designed for those
;    occasions when all the integers in arr1 will be in arr2, so these
;    indices are not needed, or when the user doesn't care about the 
;    matches from arr1.  This results in a significant increase in 
;    speed.
;
; CALLING SEQUENCE:
;    match_multi, arr1, arr2, match2, numlist=numlist, hist=hist, $
;                 reverse_indices=reverse_indices
;
; INPUTS: 
;    arr1: integer array containing unique integers
;    arr2: integer array which may contain duplicates.
;
; OPTIONAL INPUTS:
;    NONE
;
; KEYWORD PARAMETERS:
;    NONE
;       
; OUTPUTS: 
;    match2: match indices. The indices of the second list that have matches
;          in the first. 
;
; OPTIONAL OUTPUTS:
;    numlist: the number that matched each element in arr1
;    hist: the histogram of arr2, for numbers 0 to max(arr1)
;    reverse_indices: the reverse indices from the histogram function on arr2
;
; CALLED ROUTINES:
;    histogram
; 
; PROCEDURE: 
;    see code
;	
;
; REVISION HISTORY:
;    Jun-18-2002: Erin Scott Sheldon UofMich
;       Adapted from the old match2rand program
;                                      
;-                                       
;
;
;
;  Copyright (C) 2005  Erin Sheldon, NYU.  erin dot sheldon at gmail dot com
;
;    This program is free software; you can redistribute it and/or modify
;    it under the terms of the GNU General Public License as published by
;    the Free Software Foundation; either version 2 of the License, or
;    (at your option) any later version.
;
;    This program is distributed in the hope that it will be useful,
;    but WITHOUT ANY WARRANTY; without even the implied warranty of
;    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;    GNU General Public License for more details.
;
;    You should have received a copy of the GNU General Public License
;    along with this program; if not, write to the Free Software
;    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
;
;


pro match_multi, arr1, arr2, match2, numlist=numlist, hist=hist, reverse_indices=reverse_indices

  if n_params() lt 2 then begin 
      print,'-Syntax: match_multi, arr1, arr2, match2, numlist=numlist, hist=hist, reverse_indices=reverse_indices'
       print,'This works even for subsamples of subsamples'
      return
  endif 

  n1 = n_elements(arr1)

  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  ;; using histogram this way, with min=0, forces a bin for every
  ;; integer from 0 to the max. that way, we can use
  ;; arr1[i] as subscript for reverse_indices! then it
  ;; works for sub-samples
  ;;
  ;; can use lots of memory, though, if the numbers in arr1 span
  ;; a large range, since memory usage is proportional to 
  ;; sizeof(long)*span
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  hist=histogram(arr2, reverse_indices=reverse_indices, $
                 min=0, max=max(arr1) $
                )

  ptrlist = ptrarr(n1)
  numlist = lonarr(n1)
  nmatch = 0l
  for i=0l, n1-1 do begin 
      
      index1 = arr1[i]

      ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
      ;; Use of histogram with min=0 is why this works for 
      ;; subsamples of subsamples: when you do histogram
      ;; on integers with bin=1, min=0, bins are created
      ;; for each integer from 0 to max(zindex), and 
      ;; we can just subscript reverse_indices with zi. 
      ;; this could be bad for memory if all the integers
      ;; in arr1 are large
      ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

      if reverse_indices[index1] ne reverse_indices[index1+1] then begin 
          
          s2 = reverse_indices[ $
                   reverse_indices[index1]:reverse_indices[index1+1] -1 $
                   ]
          ns2 = n_elements(s2)
          ptrlist[i] = ptr_new(s2)
          numlist[i] = ns2
          nmatch = nmatch + ns2

      endif 

  endfor 

  match2 = lonarr(nmatch)
  beg=0l
  for i=0l, n1-1 do begin 
      if numlist[i] ne 0 then begin 
          match2[beg:beg+numlist[i]-1] = *ptrlist[i]
      endif 
      ptr_free, ptrlist[i]
      beg = beg + numlist[i]
  endfor 

  return
end 
